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Abstract 



r~| | Using the orbital-free density functional theory as a model theory, we present an 

Q-e analysis of the field theoretic approach to quasi-continuum method. In particular, 

r sj$ . by perturbation method and multiple scale analysis, we provide a formal justification 

for the validity of the coarse-graining of various fields, which is central to the quasi- 
continuum reduction of field theories. Further, we derive the homogenized equations 
that govern the behavior of electronic fields in regions of smooth deformations. Using 
Fourier analysis, we determine the far-field solutions for these fields in the presence of 
' local defects, and subsequently estimate cell-size effects in computed defect energies. 
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The quasi-continuum method has, in the past decade, become an important computational 
technique to study the behavior of defects in materials where a wide range of interacting 
length scales become important. The main idea behind the quasi-continuum method is a 
seamless bridging between the various length scales of interest by imposing kinematic con- 
straints on the degrees of freedom and systematically coarse-graining away from the regions 
of interest. The quasi-continuum method was first proposed in the context of empirical in- 
teratomic potentials flladmor et al.| , |1996| ), where the energy of the system was expressed 



as a non-local sum over the positions of atoms. The kinematic constraints on the posi- 
tions of atoms — degrees of freedom in the formulation — are imposed via an unstructured 
finite-element triangulation of atomic positions with full atomistic resolution in regions of 
interest, for instance at the core of a defect, and rapidly coarse-grains away to capture the 
long-range elastic effects. Apart from the kinematic constraints introduced on the degrees of 
freedom, further approximations are introduced to reduce the computational complexity of 
the formulation. The differing nature of these approximations, which include invoking the 
Cauchy-Born hypothesis in some regions of the model or introducing cluster summation rules 
in the spirit of numerical quadratures, have resulted in many different formulations of the 
quasi-continuum method. We refer to |Slicnoy ct al.| ( |1999| ); [Knap & Ortiz| (|2001|) ; |Miller fc 



[Tadmor| (|2002|) ; pliimokawa et al.j (|2004|) ; [Eidel fc Stukowski| (|2009|) and reference therein for 



a comprehensive overview of the different formulations of the method. Recent investigations 
and numerical analysis of the method ( |Shimokawa et aT| , [2004j ; |E et al.| , |2006| ; pobson ~& 
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Luskin, |2008| ; [Luskin <k Urtnci] , |2009| ; pobson et al.| , |2009| ) suggest that these approxima- 
tions can result in undesirable consequences, namely, lack of a variational structure, lack 
of stability and consistency of the approximation schemes, and uncontrolled errors in some 
cases. 

In a recent work (|Gavini et al.| > |2007a|) the quasi-continuum method was developed for elec- 
tronic structure calculations using orbital-free density functional theory (OFDFT). OFDFT, 
which is an approximation to the widely used Kohn-Sham formulation of density functional 

1964| ; [Kohn fc Sham| , |1965| ), describes the ground-state energy 



theory (Hohcnbcrg fc Kohn 



of the system as an explicit functional of electron-density and is valid in materials systems 
whose electronic structure is close to a free electron gas (cf. |Parr fc Yang| ( p. | ) ; [Wang 



fc Teterj ([1992]); |Smargiassi fc Madden| (|i"994|) ; |Wang et alj fll99l| [T999|) for a comprehen- 
sive overview). The quasi-continuum reduction of OFDFT was achieved using a real-space 
local variational formulation, and a coarse-graining of the electronic-fields and positions of 
atoms — degrees of freedom in the formulation — through kinematic constraints imposed using 
nested finite-element triangulations. An important difference in the mathematical structure 
of quasi-continuum formulation for OFDFT in comparison to empirical interatomic potentials 
is that OFDFT is a local field theory as opposed to the non-local description of extended 
interactions in empirical potentials. A local field formulation, as in the case of OFDFT, 
admits quadrature approximations to further reduce the computational complexity with- 
out introducing the undesirable consequences characteristic of conventional quasi-continuum 
formulations. 

In the prequel to this article (flyer fc Gavini| , |2010|) , we suggest a field formulation for 
commonly used interatomic potentials, where the extended interactions in these potentials 
are reformulated into a local form by constructing partial differential equations (PDE's) 
whose Green's functions correspond to the kernels of the non-local interactions. We further 
demonstrate that the quasi-continuum reduction of these field formulations is variational, a 
consistent numerical approximation, and provides significantly better accuracy than previous 
formulations. Moreover, the field formulation of interatomic potentials provides a unified 
framework where the quasi-continuum reduction is solely a numerical approximation scheme 
irrespective of the field theory used to describe the system — density functional theory or field 
theories that represent interatomic interactions. 

In the quasi-continuum reduction of field theories ( [Iyer fc Gavini| , [2010| ; |Gavini et a 



2007a| ) , the various fields that appear in the formulation are decomposed into predictor fields 



and corrector fields. The predictor fields are computed by performing a periodic calculation 
using the Cauchy-Born hypothesis, and the corrector fields are subsequently computed from 
the variational formulation. For smooth deformations which do not depend on the atomic- 
scale, plane et al.| (|2002|) show that the various fields are given, to the leading order, by 
a periodic calculation using the Cauchy-Born hypothesis. Hence, in regions away from the 
defect-core it is expected that the predictor fields are good approximations to the fields. 
Thus, the corrector fields are represented on a finite-element triangulation which is subatomic 
near the core and coarse-grains away to become superatomic, and this constitutes the quasi- 
continuum reduction of field formulations. 

The representation of the corrector fields on a coarse-grained triangulation is valid un- 
der the hypothesis that corrector fields do not exhibit oscillations on the atomic-scale. In 
this work, we provide a formal justification for this hypothesis. We conduct our analysis 
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in the framework of OFDFT and latter comment on other field theories. We first use the 
perturbation method to find the governing equations for the corrector fields corresponding 
to a weak defect. While the defect plays the role of a source (forcing function), the coeffi- 
cients of these governing equations are given by the unperturbed (predictor) electronic fields 
and hence oscillate on the atomic-scale. For homogeneous deformation, the unperturbed 
electronic fields are given by periodic calculations with respect to the unit cell at atomistic 
scale; for a smooth macroscopic deformation, the unperturbed electronic fields are generally 
unknown. Motivated by the thermodynamic limit ( Blanc et al. , 2002 ; |Garcia-Ccrvcra et al 



2007|) , we nevertheless hypothesize the unperturbed fields are given by periodic calculations 
with respect to the local atomistic lattice. Further, since the unperturbed electronic fields 
oscillate at the atomistic scale which is much smaller than the macroscopic supercell and 



the considered macroscopic perturbation, we employ the multi-scale analysis ( |Cioranescu &: 



Donato , | ) to find the corrector electronic fields. In particular, we demonstrate that the 



corrector electronic fields to their leading order and first order are independent of the lattice 
parameter, and hence do not exhibit atomic-scale oscillations. However, this result shall be 
interpreted with caution near the defect-core since in reality a defect, e.g., a vacancy or an 
interstitial, is localized at the atomic-scale. Further, we derive the homogenized equations 
for the corrector fields. These homogenized equations turn out to be a second-order linear 
system of PDE's. By Fourier method, we find their Green's functions explicitly, which show 
that the correctors fields in OFDFT, corresponding to electrostatic potential and electron 
density, decay exponentially. Additionaly, we compute their solutions for a situation repre- 
sentative of a vacancy in an infinite crystal. Using these solutions, we analyze the cell-size 
effects arising from a computation on a finite domain and estimate the domain size required 
for achieving chemical accuracy in vacancy formation energy. Our results show that a cell-size 
of the order of 1,000 atoms is required to attain a converged value for the vacancy formation 
energy in aluminum, which is much larger than the cell-sizes that are commonly used in 
numerical simulations. This estimate is in agreement with a recent cell-size study in |Gavini 



et al.| ( |2007a| ). We further note that the cell-size effects are likely to be more significant for 



defects like dislocations where the decay in elastic fields is much slower. 

The remainder of this paper is organized as follows. In section [2| we formulate OFDFT 
with Thomas-Fermi- Weizsacker kinetic energy functionals, and present the problem defini- 
tion and the assumptions made in this analysis. In section |3| we discuss the perturbation 
analysis of the corrector fields, and present the multi-scale analysis of these fields and derive 
the homogenized equations in section [|. In sections |5| and ^ we derive the Green's functions 
of the homogenized equations and compute their solutions for a spherical defect representing 
a vacancy. In section ^ we comment on the extension of this analysis to other flavors of 
OFDFT which use non-local kernel energies and other field formulations representing empir- 
ical interatomic potentials. We finally conclude in section |] providing an outlook. 

2 Problem definition 

Consider an infinite crystal with lattice points given by C v = rjC and C = {X^i=i '■ 
vi,V2,v 3 G 2Z}, where r] « 1 denotes the fixed lattice parameter and ei,e 2 ,e 3 e IR 3 are 
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the rescaled lattice vectors satisfying e 3 ■ (ei x e 2 ) = 1. We refer to 



U v = rjUo, 



Un 



i=i 



1 !i 



as the unit cell and the rescaled unit cell, respectively. Let Y = (—1, l) 3 be a macroscopic 
supercell such that it overlaps with an integer number of unit cells, Z be the charge at each 
nucleus measured in units of electron charge, and y : Yq — > Y be a smooth macroscopic 
deformation that carries a reference point xo 6 Yq to a new point y(xo) £ Y. In this 
work we are interested in macroscopic deformations that are independent of rj. We assume 
that the nuclei follow the Cauchy-Born rule, and hence the nuclear charges in the deformed 
configuration are given by 



y( x o)), 



where 5 is a regularization of the Dirac distribution that represents a unit nuclear charge. 

To present our ideas we consider the energy of a system described by OFDFT. We remark 
that the ideas presented here are general and can be equally applied to any field theory, for 



instance, fields theories that describe empirical interatomic potentials discussed in [Iyer & 



(|2010|) . In density functional theory, the energy of a material system is given by 



E(u,b y ) = T s (u) + E x 



u 



E H (u) + E ext (u, b y ) + E zz (b y 



(1) 



where u denotes the square-root electron density, T s denotes the kinetic energy of non- 
interacting electrons, E xc denotes the exchange and correlation energies that account for 
the quantum mechanical effects, and Eh, E ext , E zz denote classical electrostatic interaction 
energies between electrons and nuclei. In OFDFT, T s is approximated by explicit functional 
forms of electron density as opposed to the Kohn-Sham approach where it is computed 
exactly within the mean field approximation. A simple choice for this approximation is the 
Thomas- Fermi- Weizsacker (TFW) family of kinetic energy functionals ( Parr fc Yangj , |1989| ): 



T s (u) 



C F j 



(2) 



Y 



where < A < 1 is a parameter and Cp = ^(37r 2 ) 2 ' 3 . More accurate kinetic energy 
functionals have been proposed in the past decade which account for the linear response of a 
uniform electron gas. For clarity we postpone our analysis of these functionals to section [7|. 
By choosing the Thomas-Fermi- Weizsacker functionals fl2|) for kinetic energy, and following 
the real-space formulation of OFDFT proposed in |Gavini et al.| ( |2007b|) , we express the total 
energy of the system as 



E(u 



, 0; 6 y ) 



f{u) + ±\Vu\*- l -\V<P\ 2 + {u* + b y ) 



<ix, 



(3) 



where f(u) = Cfu 10 ^ 3 and (f> denotes the trial electrostatic potential. In the above, we ignore 
exchange and correlation energies and comment on them in section 0. The ground state of 
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(0, u, y) is determined by the following min-max problem 



£tot(0) := min \S(y) := min max E(u,4>;b y )\, (4) 

where 

3^ : = {y G x : Vy periodic in F }; 
U(by): = {u E H l per (Y) : ^_(n 2 + 6 y )dx = 0, u > 0}. (5) 

In the above definitions, x is a suitable function space that admits minimizers of £(y). In this 
analysis, since our focus is to derive and analyze the far-field behavior of the displacement 
and electronic fields, we restrict our attention to a local minimizer of £(y) in y. 

Let (0y,w y ) be a solution of the min-max problem for a smooth deformation y 6 ^ 

£(y) = E(uy, y ; b y ) = mm max E(u, 0; b y ). (6) 

The existence of a solution for the saddle point problem (^|) can be established following the 
ideas in |Gavini et al.| ( [2007b|) , where the analysis was carried out in a non-periodic setting 



with Dirichlet boundary conditions on a bounded domain. We remark that the arguments in 
Gavini et al.| ( 2007b|) can be appropriately modified to the periodic setting, and these details 



are not discussed in this article to maintain our focus on multi-scale analysis. We also refer 
to |Licb| ( |1981|) for results on the existence and uniqueness of solutions for various flavors of 



OFDFT. 

It is clear from the definition (|3|), if (0 y ,u y ) is a solution to the min-max problems in 
equation (||), so is (0 y + c, u y ) for any c G IR. By the standard first- variation calculations 
it follows that there exists a solution to the min-max problem in equation (^[), denoted by 
(0 y , u y ), satisfying 

{A0 y + (u 2 y + b y ) = on F, 

-AAu y + /'(%) + 2%0 y = on Y, (7) 

subject to: u y G U{b y )] <f) y G Hp er (Y). 

Note that, in the above equation, the Lagrangian multiplier associated with the constraint 
in equation (|) 2 has been absorbed into the electrostatic potential y . Thus, the solution </> y 
to problem ([7|) no longer allows an arbitrary additive constant (cf. |Catto et al.| ( |l r ] ) for 
further discussion on this point). 

We now discuss the nature of the solution (<fr y ,u y ) to problem (|7]). First we assume 
a homogeneous deformation with Vy = F G IR 3x3 on Y . Consider problem (|7|) on the 
deformed unit cell Fq^ 

( A0 + (u 2 + b y ) = on F U V , 

< -AAw + f'(u) + 2u(f) = on F ^, (8) 

[subject to: j ¥ ^ v ^{u 2 + 6 y )rfx = 0, u, G Hp er (F U v ). 
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For f(u) = ^,(3-7t 2 ) 2 / 3 m 10//3 , |Catto et al.| ( |1998| ) have shown that the periodic extension of the 



solution to problem @ with respect to the period FqU v , denoted by ((f)*, u*), is the solution 
to problem (0): 



x 



w y (x)=w*(x) VxGF Fo- (9) 



For future convenience, we denote by 

p (F o ,x) = 0*(r/x), M P (F ,x) = u*{rpc), (10) 

where the subscript p signifies that x H- (0 P (F O , x), u p (F , x)) are periodic with period equal 
to the rescaled unit cell F q Uq. It is worthwhile noticing that (4> p ,u p ) are considered as being 



defined by the exact solutions to the unit cell problem (||) through equation flTo"D , instead of 
the solutions in the thermodynamic limit discussed in [Blanc et al.| ( [2002|) . By equations ([|) 
and (IToT) , we have 



'y 



x) = 0p(F o ,x), tt y (x) = tt p (F ,x), (11) 



where x = x/77 denotes the fast variable in the subsequent homogenization calculation. We 
remark that, in Section 4, equations ([LTD and the fact that 77 << 1 compared with the macro- 
scopic supercell Y will be used to derive the homogenized equations for the corrector fields. 
This homogenization limit is not the thermodynamic limit where the nuclei are assumed to 
locate at r/C fl Y and 77—7-0. Trying to couple the homogenization limit and the thermody- 
namic limit encounters difficulties, and we will not address this issue in this paper. Further, 
we identity the solution to problem (Mj is also a solution of the min-max problem 



/r A 1 

/(«) + o I Vm I 2 - o l V< ^! 2 + + 6 y)0 
f u„ L 2 2 



dx (12) 



subject to the same constraints as in equation (|8]). Here and subsequently, ) = ^^~7q) Jq( 
denotes the averaged value of the integrand over the domain Q. In terms of the solutions 
(<f) p , u p ) to the unit cell problem, we define the following quantities for future use 

a(F ) = f Wp(F ,x)dx, /3(F ) = f p (F o ,x)rfx, 

7 (F )= / [i/"K(F ,x))]rfx + /3(F ). (13) 

7 F C/o 2 

As in classical continuum mechanics, all the functions W, a, /3,j : iR 3x3 — > IR satisfy the 
material frame indifference and material symmetries 

W(RF ) = W(F ) VR G SO(3) & F G IR 3x3 with detF > 0, 

W(F U) = W(F ) VH G ^(F ), (14) 

where 50(3) consists of all rigid rotation matrices, G(F ) is the point group associated with 
the Bravais lattice J-qC, and W can be replaced by a, (3, or 7 in equation ([141) . Equation ( |T~4|) 
can be verified directly from the definitions (^2|) and (|T3|). 
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We now consider the case when the deformation y has a smooth macroscopic deformation 
gradient F : Y — > IR 3x3 on the current configuration 

F(x) = F(y" 1 (x)) V x e Y, F(x ) = Vy(x ) V x e Y . (15) 

A priori, for this case, we have no knowledge on the solution to (0). Since rj « 1, motivated 
by [Blanc ct~aT| ( [2002| ) we hypothesize that the solution to (|7|) is given by 

y (x) = p (F(x),x), M y (x) = w p (F(x),x), (16) 

and the elastic energy is given by 

£(y) = E(u y , y ; by) = J W{F{x))dx, (17) 

where (<f> p , u p ) are defined by the exact solutions to the unit cell problem (^) through equa- 
tion (0), and W : IR 3x3 -> IR is the elastic energy density on the deformed configuration Y 



given by equation (|12|) 



The solution to the outer minimization problem (|]) may not be unique, and throughout 
this work we will restrict our attention to local minimizers that satisfy the Euler-Lagrange 
equation corresponding to the energy in equation ([17]), which is the familiar equilibrium 
equation of elasticity 

DivS(Vy*) = on Y , (18) 

where 

dJW 

[S(Vy*)] pi = ^r-(Vy*), m = det(F). 



<9S 

Note that in equation ([Taj, S is the first Piola-Kirchhoff stress and DivS 
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3 Perturbation analysis 

We now consider the effect of defects on electronic fields ((f), u). A defect breaks the lattice 
symmetry which in effect is a perturbation of the forcing term, b y , in equation ([?]). Thus, 
we replace the forcing term in equation (0) by a small perturbation of b y : b y = b y + eb c 
with e << 1 and consider b c to be independent of the lattice parameter 77 which allows us to 
subsequently pass to the homogenization limit in section |J If b c has a compact support, this 
perturbation can be interpreted as a weak local defect, formed by slowly reducing the charges 
on the nuclei in a macroscopic region, in an otherwise perfect crystal undergoing a smooth 
deformation. We are interested in calculating the influence of this perturbation (defect) on 
the ground state of OFDFT and, in particular, on the total energy. As in equation @, the 
ground state of the system is governed by 

£tot(b c ) :=mm\S e (y;b c ) := min max E(u, 4>;b y + eb c )\ . (19) 
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Note that if b c = 0, i.e., the system is unperturbed, then £tot(b c = 0) is equal to £f O t(0) in 
equation (||). 

We solve the above problem approximately by perturbation method. We first consider 
the inner min-max problem in equation (|l^) for given y G y. Let 

<P £ = <P y + e<p c G Hp er (Y), u £ = u y + eu c G U(b y ), (20) 

be the solutions, where (0 y , u y ), the solutions to the unperturbed problem (^), are referred to 
as the predictor fields in the quasi-continuum formulation, and (0 C , u c ) are referred to as the 
corrector fields (|Gavini et al.| , |2007a| ). Inserting equation ( p0[) into equation (^) 2 , we obtain 



the charge neutrality constraint 

J [2u y u c + eu\ + 6 c ]dx = 0. (21) 
Inserting equation (^) into equation (^), we expand the energy as 

E(u £ , 4> £ ; by + eb c ) = E(u y , <p y ; b y ) + e J | [f'(u y ) - AAu y + 2u y y ] u c 

+ [A(f>y + u 2 y + b y ]4> c + 6 c y |dx + £ 2 E 2 (u c , C , y; b e ) + o(e 2 ) 
= E(u y , y ; by) + e J 6 c y dx + £ 2 E 2 (u c , C , y; b c ) + o(e 2 ), (22) 

where the second equality follows from the Euler-Lagrange equations in (|7|) for (<f) y , u y ), and 

/r 1 a A 1 i 

—f"{uy)u c 2 + -|Vw c | 2 - -|V0 C | 2 + (2u y u c + 6 C )0 C + (f) y u 2 c (ix. (23) 

Neglecting o(e 2 )-terms in equation (|22|), by the inner min-max problem in equation (19) we 
arrive at the following min-max problem for (u c ,<f) c ): 

£2(y;b c ) := mmmaxE 2 (u c ,(f) c ,y;b c ) (24) 

U c (j) c 



subject to the constraints (cf. equation (|2lD) 

C G H^Y), u c eH x per {Y), J^2uyu c + b c )d* = 0. (25) 



We remark that the zeroth and first order terms in equation ([22] ) are absent in the min- 
max problem fl24T) since they are independent of (<p c ,u c ). By the standard first- variation 
calculations, we show that the Euler-Lagrange equations for (<fr c , u c ) associated with the 



min-max problem (pi) are 

A0 C + (2u y u c + b c ) = on Y, 

■\Au c + (/ (rt y ) + 2(py)u c + 2uy(f> c = on Y, 

where, as in equation (0), we have absorbed into the potential C the Lagrangian multiplier 
associated with the last constraint in equation (p5|), which is a constant independent of x. 



We further notice that equations (^6|) can be obtained by linearizing equation ([F]) near the 
solutions (0 y ,M y ). We remark that although the perturbation analysis was conducted under 
the assumption of weak local defects, the perturbation expansion given by equation (^) is 
a reasonable assumption in regions away from defects that are not necessarily weak. This 
follows as the perturbations in electronic fields decay away from the defect core due to the 
elliptic nature of the PDE's, and the governing equations for corrector fields will subsequently 
be valid in these regions. 



4 Homogenization 



We now turn towards establishing certain properties of the corrector fields which play a fun- 
damental role in the construction of quasi-continuum reduction of field formulations proposed 
in pavini et al.| ( |2007a| ); |lycr fc Gavini| ( |2010| ), and provide a formal mathematical justifica- 
tion for the method. Before proceeding to details, we notice the following useful identity. Let 
/(x, x) be a smooth function which is periodic in the second variable x with period F(x)[/o- 
If 7] « 1, we have the identity (cf. e. g. |Cioranescu &: Donato| ( |1999| ), Chapter 2), 



/ /(x, — )<ix = / +- /(x, x)<ix<ix + o(l). 

JY V JY J F(x)C/ 



(27) 



Since the unperturbed solutions (0 y ,w y ) given by equation (Q) oscillate at the atomic 
scale-?], presumably the corrector field solutions (<f) c , u c ) to the governing equations in ( f2~6|) 
oscillate at the 77-scale as well. In this section, we determine the order of this 77-scale oscillation 
in the corrector fields (<fr c , u c ) and whether this atomic-scale oscillation is important to the 
leading order in energy. Further, we determine the homogenized equations that govern the 
macroscopic behavior of these corrector fields. To this end, following the method of the 
multiple scale expansions, we assume 

. c (x) = $(x,x) + 77#(x,x) + --- , 
(x) = u°(x,x) +?7uJ(x,x) H , 

where x = x/7/ is the fast variable, 0*(x, x), u*(x, x) {i = 0,1,- 
periodic in the fast variable x with period F(x.)Uq. Replacing (0 y , u y 
of equation (|l^), we rewrite E2 in equation (|23|) as 



(28) 



■ ) are assumed to be 
by the right hand side 



E 2 (u c ,(fi c ,y;b c 



, 1 r<> 

k "2 



4>p)u c 2 + ^|Vw c | 2 - ^|V0 C | 2 + (2u p u c + b c 



dx. (29) 
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Inserting the multiple scale expansion ( |28|) into equation (|29|), we have 



E 2 {u c ,<j) c ,y;b c 



1,1. 



~|-V^° + V x 0° + V x 0*| 2 + (2 MpM ° + b c )<P° c rfx + 0(r/) 
z 77 J 



2r7 2 



r 



A|V X «°| 2 - IV* 



/ 0|2 



y^rc l 



1 

V Jy 



+ 



AV X «° • (V xM ° + V*^) - V x 0° • (V x 0° + V*<#) 
(\f"(u p ) + (f) P )(u° c ) 2 + ^|V x w° + VscuH 2 
-^|V X 0° + V x 0*| 2 + (2u p u° + 6 C )0°] rfx + 0(77). 



dx 



(30) 



We neglect the higher order terms of 0(r]) in equation (|30|) . Since 77 << 1, we consider the 
min-max problem (24) first for the leading ^ -terms in equation (|30|), which is given by 



mm max 

u° <f>° 



A| V x u 



1 2 



IV*. 



t0|2 



c?x. 



It is clear that a solution to the above problem necessarily satisfies 
V x w°(x, x) = and V x 0°(x, x) = 0, 



(31) 



which means that (jP c (x, x) and u° c (x, x) are independent of the fast variable x and hence can 
be rewritten as 



u°Jx) 



and 



x . 



(32) 



This shows that the leading order terms in the corrector fields do not exhibit atomic-scale 
oscillations, and thus the corrector fields can be resolved accurately on length scales larger 
than the lattice parameter. This key result formally justifies the coarse-graining of corrector 
fields introduced in the quasi-continuum reduction of field theories. 

Further, in account of equation (|3~TD , the ^-terms on the right hand side of equation ( Bof ) 
vanish. Finally, we consider the ^--order terms on the right hand side of equation ([30]) which 
represent the leading order terms in the multiple scale expansion of E 2 . Using equation ( E7p 
we can rewrite equation (BO) as 



E 2 (u c ,4> c ,y; b c 




A 



V x n 



1 2 



i0 1 2 



+ 



'Y J F(x)C/ 

If- . 

'Y J F(x)t/ 



cZxdx 



-/"K(F(x),x)) + p (F(x),x) )\u 



,0 1 2 



+27i p (F(x),x) M °0O + 6 c 0O 



dxdx 



+ 




YJ F(x)C/ 



^(2V X «° + V x <4) • V x ^ - (2V X 0° + V x <^) 



(33) 



dxdx.. 



10 



Since mJ(x, x) and 4>l(x, x) are periodic on F(x)£/" for every x, from equation ( |3"2"D we have 



V x m° • V x M^rfx = 



and 



F(x)t/ 



/ Vx0°-Vs^dx=O VxgF. (34) 

JF(k)Uo 



From the min-max problem (|24"D, we maximize the expression in equation (Q) over admissible 



!>* and minimize it over admissible ul, and obtain 



V x ^(x,x) = and V x 0^(x,x) = O. (35) 
Thus, it follows that the corrector fields do not exhibit atomic-scale oscillations up to the 



second order terms in the multiple scale expansion (28). Further, from equation (|35|) , the 
last term on the right hand side of equation ( |3~3|) vanishes, and by equations ([13|) and ( p2[) 
we identify the first two terms in equation (|33D as 



y 



^IVx^I 2 



/ 1 2 



7 (F(x))|^| 2 + 2a(F(x))«»^ + rfx =: y; 6 C ),(36) 



,0i0 



0/„,0 lO 



where the 7 and a are defined in equation (|13|) . In conclusion, from the min-max problem 
in equation (^4|), assuming the multiple scale expansion given by equation (|28|), and keeping 
only the leading order terms, we have 



subject to 



£2(y;b c ) « minmax£ 2 (w c ,0 c ,y;6 c 

u° 4>° 



(37) 



(3? 



where the constraint on u° c follows from equations (|25|) , (13), and neglecting higher order 
terms in equation ( p8|) . Equations (j37|)-(|3"8|) constitute the governing equations for the cor- 
rector fields in their leading order. 

We now proceed to derive the governing equations that describe the elastic response of 
the defect. From equations (|22[) and (|37|) we see that £ e : y — >• IR defined in equation ( |i~9D 
is given by 



£ £ (y;b c ) = E(u y , y , by) + e J 6 c y rfx + e 2 £ 2 (y; b c ) + o(e 2 ) 
W(F(x))dx + e f b c (f) y dx + e 2 £ 2 (y;b c ) + o(e 2 ) 



y 



(39) 



where in the second equality we have used equation ([17]) for E(u y , <p y , b y ). Let y* : Y — > Y 
be the unperturbed minimizer of the outer minimization problem in ([|), 

F*(x ) = V X() y*(x ) Vx eF , F*(x) = F*(y*- 1 (x)) Vx G Y, 

and, parallel to equation (|2H|), let 

y £ = y* + ^y c (40) 
11 



with y e G y being the minimizer of the outer minimization problem in (|19l). As we are 
interested in the elastic fields created in response to the perturbation b c , and not the con- 
figurational force on b c , we hold the pull back of b c on to the reference configuration fixed. 
To this end, we define b c (-x ) = 6 c (y*(x )) Vx G Y as the pull back before introducing 
the perturbation in the deformation field. Subsequently, for any infinitesimal perturba- 
tion of the deformation given by equation (|40|), b c in the current configuration is given by 



6 c (x) = b c (y £ *(x)) Vx G Y. Further, let 

F £ (x ) = V X() y £ (x ) VxqGFo, F £ (x) = F £ (y £ "(x)) Vx G Y, 

and 

d 2 TW 8 713 

[C(V xo y*)U- := ^-^r-(V xo y*), [B(V xo y*)] pi := ^(V^y*). (41) 

Since 

/ W(F £ (x))dx = / J(F £ (x ))W(F £ (x ))dx , 

JY JYq 

we have 

/ PF(F £ (x))dx = / J(V Xo y*)^(V Xo y*)dx + £ / V Xo y c ■ S(V Xo y*)dx 
Jy Jy Jy 

+e 2 I ^V Xo y c • C(V Xo y*) • V Xo y c rfx + o(e 2 ). (42) 

JYo 

Further, by equations and (p~3| ) we have 

/ b^dx « / 6 c (x)/ p (F £ (x),x)rfxrfx= / 6 c (x )J(F £ )/3(F £ )dx 
Jy Jy J f^(x)u Jy 

= [ 6 c (x )J(Vy*)/3(Vy*)dx + £ / 6 c (x )V X() y c ■ B(Vy*)dx + o(e) 

Jy Jy 

« / fe^dx + e / 6 c V X() y c • B(Vy*)rfx . (43) 

Replacing y in equation fl39f) by y £ given by equation (f4(|), expanding and keeping terms up 
to e 2 , by equations (f4"2|) and (|43|) we obtain 



£ £ (y* + ey c ; & c ) « / J(V xo y*)W(V xo y*)dx + e [ V xo y c • S(V X() y*)dx + e [ b c <f> y *dx 
Jy Jy Jy 

+e 2 J {^V Xo y c • C(V Xo y*) ■ V Xu y c + 6 c V Xo y c • B(V Xo y*)}dx 

+e 2 S 2 (y*;b c ) + o(e 2 ), (44) 

Since y* is a local minimizer satisfying equation ([18]), it follows that J Y V xo y c -S(V xo y*)<ixo = 
0. We further neglect o(e 2 )-term in equation (|TT|). Finally, the outer minimization problem 
given by equation (|T9D reduces to a minimization problem on y c and is given by 

S el (b c ) := min / {4v xo y c ■ C(V xo y*)V xo y c + 6 c V xo y c • B(V xo y*))dx . (45) 
y c ey J Yq < 2 J 
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Thus, a minimizer y c satisfies the following Euler-Lagrange equation which constitutes the 
governing equation for the elastic response in the presence of a defect 



Div[C(V Xo y*) V Xo y c + 6 c B(V Xo y*)] = on Y . (46) 

An important quantity in the study of defects is the defect formation energy, which is 
defined as the excess energy in the system with a defect measured from a reference state of a 
perfect crystal consisting of same number of particles — in this case the number of electrons 
and nuclei. In the framework of the present study, it is given by 



£d(b c ) := [£tot(b c ) - £tot(0) - e J 6 c y .dx]/e 2 . 

From the previous discussions, the defect formation energy (defect energy) can be expressed, 
to the leading order, as the following min-min-max (saddle point) problem: 

8 d (b c ) « minminmax f (4v xo y c ■ C(V Xo y*)V xo y c + 6 c V xo y c • B(V Xo y*))dx 
y c u o 4,0 J Yo I 2 ) 

+ / {7(F*K| 2 + ^|V x u°| 2 - i|V x ^| 2 + 2a(F»° + 6 c 0°}dx (47) 

subject to 

fieH^Y), uleHlXY), J(2a(^)u° e + b e )dx = t y c ey. (48) 

Associated with the above min-min-max problem, the Euler-Lagrange equations are the 
elasticity equation ([46]) for y c on the reference configuration Y and 

A0O + 2a(F*)u° + b c = on Y, 

-AAu° + 2 7 (F*)u° + 2a(F*)<#> = on Y 



for (0°,m°) on the current configuration Y. Note that the elasticity problem ([46]) for y c is 
not coupled with the equations for (0°, u° c ). In terms of the solutions (y c , 0°, u° c ) to equations 
(^6|) and (p49|) , the defect energy can be written as 

£d(b c ) w 4 / ^ V y c • B ( v y*) + 4 / (50) 



JYo z JY 



5 Far fields 

In this section we determine the far- field behavior of the fields (0°, u° c , y c ) from the governing 
equations in (|49|) that will aid in determining the optimal coarse-graining rates for these 
fields. In this analysis, we assume b c is continuous, bounded and supported within the ball 
B ro = {x : |x| < r }. Although the analysis in the previous section was performed on 
the supercells Yq and Y, we note that the results of the analysis are independent of the 
supercells and thus to determine the asymptotic behavior of the corrector fields we assume 
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Yq = Y = IR 3 . We first calculate the far field behavior of for a homogeneous 

deformation with V xo y* = F G IR 3x3 on IR 3 . In this case, F* = F on IR 3 as well; 
C(V Xo y*), B(V xo y*), a(F*), 7(F*) are constants on IR 3 and we drop their dependence on 
Fo in notation. Further, the periodic boundary conditions in equation ([48]) shall be replaced 
by appropriate decay conditions at the infinity. Dropping the subscript c and superscript 
in ((j) c ,u° c ,y c ) in equations (|49|) and (fl6|), we rewrite our problem for (</>, u,y) as 



A<p + 2au + b c = on IR 3 , 

-A Ait + 2 7 u + 2a(f) = on IR 3 , 

Div[CV Xo y + 6 c (F x )B] = on IR 3 , 



(51) 



subject to 



|n(x)|, |y| —7-0 as |x| — > +oo, 



(2cm + b c 



(52) 



We now address the solutions of the first two of equation (^T|). Taking Laplacian of 
equation (|5l"]) 9 and inserting into equation (151]),, we obtain 



2 1 
AAu — -j An + = 6 



where Z > 0, Re(Zi) > 0, 



7 



Z? A' Z 



4^ 
A 



> 0, 6 



on iR 3 



2ab r 



(53) 



A 



(54) 



The constants l\, Iq determine the asymptotic behavior of the fundamental solution at 
the infinity. Since equation (Q) is linear, we express its solution as 



u(x) = (£ u *6)(x) = / E u (x-x!)b(-x!)dx! 

Jr 3 



(55) 



where E u is the fundamental solution satisfying 



(AA-^A + -j)^ = *(<)), 

and 5(0) is the Dirac distribution. We find the fundamental solution E u by Fourier analysis. 
Solving the algebraic equation 



4 2 „ 1 

X + 12 X + 74 = °. 



(56) 



we obtain two roots k± with Im(K±) > and satisfying 



1 + 1 1 



(57) 
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The two other roots with Im(K-i-) < are discarded as they will correspond to exponentially 
growing solutions in (u,</>), defined subsequently, and do not satisfy the decay conditions 
imposed in (|5*2"|). By Fourier analysis, we have 



EM 



(2tt)3 J r3 (4 - K 2 _) 



\ 2 -k 2 



I 2 - kI 



exp(zk • x)dk. (5? 



We are therefore motivated to consider the fundament solution of the operator 

(A + k 2 )G = 5(0) (59) 

for some with Im(ft) > 0. By the standard method (cf. |J acksonj (|1999| ) page 243), we 

have 



exp(ire|x|) 



G(x, k) 



if Im(«) ^ 0, 



_ A exp(^|x|) _ g cxp(^|x|) ifIm ( K ) = 0) 



47r|x| 

4-7r|x| ^ 4-7r|x| 

where A + B = 1. In Fourier space, equation (^9|) can be rewritten as 



(60) 



G(x, k) 



-1 



(2tt) 3 7 r3 |k| 2 - K 1 
If j^l 7^ / 0; i- e - ; K + 7^ from equations (|58l)-(|6l"1) we have 

-1 



exp(zk • x)dk. 



EJx) 



(k 2 - K 2 



G(x, k + ) — G(x, K. 



(61) 



(62) 



If \ — l > 0, i.e., 4 = K - = ~~ ^) sending k + to k_ in equation (|62|) we obtain 



EJx) 



if 7 > 0, k — i/l 



^d K G(x,*) = t e M-\x\/lo) 
^ K G(x, K ) = ^exp(z|x|// ) + ^exp(-z|x|// ) if 7 < 0, k = l/l . 



(63) 



Further, by the second of equation (0), the associated potential is given by 

/ £ (x - x>(x')(Ix', 
Jr 3 



(x) = — Ah — — u 

2a a 



where 



^(x) = ^[AK(x)-|k(x)]. 



(64) 



(65) 



We remark that the above formal calculations can be rigorously justified (cf. e. g. |Rudin 
( |19911) , chapter 7). 

Note that the last of equation (|52"|) requires 



( F «-6) = 0. 



(66) 
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If J 6(x)dx 7^ and 6(x) has a compact support in a ball around the origin, then for large |x| 
the solution u(x) is well approximated by the Green's function E u in equation fl62f), which 
is not integrable if rz± are real numbers as J R z cxp j^ x ^ <ix is not integrable for real fc. We 
therefore conclude that k± should be both nonreal numbers. This is possible for the following 
three cases: 



1. 7 > and Iq > l±. In this case, all roots of equation ( p6|) are pure imaginary. By 
equations fl62|) and we have 



K(x) 
EJx) 



Att{k\ - K 2 _)\x\ - 

1 

- «i)|x| 



exp(m + |x|) — exp(m_|x|) 
C + exp(z/t + |x|) — C_ exp(m_|x|) 



where 



/~y ^ / 2 2 A 2 

C+ = —(-K ± - 72 ) = — 



2a 



if 2a 



(67) 



(6* 



2. 7 > and Zq = This is the first case in equation (|63|) and we have 



K,(xl 



/r 



57T 



exp(-|x|// 



-A ,1 



167ra Iq r 



exp(-|x|// ) 



(69) 



3. /o < |/i|. In this case, all roots of equation (|~~|) are nonreal and the fundamental 
solutions are given by (|67| ) as well. 



To verify the constraint (|6q) , we integrate equation ( p3|) on the ball with radius N, and 
by the divergence theorem arrive at 



/ (VAu - -jVu) • MS + / (-jM-6)dk=0, 



(70) 



where x = x/|x|. Sending iV — > +oo, we arrive at equation ( |66|) since the first term in the 
above equation vanishes for expressions in equations (|67j ) or (|69]). 

Finally, we remark that the solution to the last of equation (|5lD is given by the classic 
theory of linear elasticity (cf. e. g. |l\lura | ( |1987| ), chapter 1). 

Equations (|55D, (|64l), (|67j ) -(|69"l ) and the theory of elasticity determine the far-field behavior 
of (u,<f),y), where the perturbation b c plays the role of a source. For a continuous bounded 
b c supported within the ball B ro = {x : |x| < ro}, we have that for some R > and some 
C, k > 0, 



l u ( x )| < Cexp(— k|x|), 



(x)|<Cexp(-K|x|), |y(x )|< 



Vlxl > i2. 



(71) 



With the above estimates on the far-fields, we continue our solutions to equations fl5l|)-(|5! 
for a particular example in the next section. 
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6 Defect energy and cell-size effects 



In this section we study how the defect energy depends on the size of the supercell. For 
simplicity, we assume that the supercell is the ball -Br = {x : |x| < Ro}, the coefficients 
l > li > and thus both the roots k± in equation ([57]) are pure imaginary numbers. We 
denote by 



k 




7 =F a/7 2 - 4a 2 A 



A 



> 0. 



(72) 



Below we solve equations (^T|) for the corrector fields (u, 0, y) with 

if |x| < r , 
if Ixl > r , 



6(x) 



where p G IR is a constant, r < Ro describes the length scale of the defect representative of 
a vacancy We apply the Dirichlet boundary condition 



u x 



0. 



A , . 7 , . 
— Au x - -u(x) 
2a a 



y = on &B Bo , 



(73) 



where <j G iR is a constant determined by the constraint (|S6|). 

We first consider the electrostatic contribution of the defect energy, i.e., the second term 
on the r.h.s. of equation fl50l). By symmetry, we have u = u(r) with r = |x|. Therefore, 
equation (|53|) can be rewritten as 



d 4 



2 d 2 



dr 4 ™ > 2 



I 2 dr 2 



-ru 



ru 



rb 



V0 < |x| < R . 



From the theory of ordinary differential equation, we obtain 

' prl\ + Ciexp(k + r) + C 2 exp(/c_r) 

+C 3 exp(-k + r) 
C5exp(fc + r) + C*6exp(/c_r) 

+C 7 exp(-k + r) 



ru\r) 



C4exp(— fc_r) 
C 8 exp(— k_r) 



if r < ro, 



if r > r , 



where the constants C, (i = 1, ■ • • ,8) are determined by the analyticity of ^(x) at r = 
(which implies u(r) is an even function, i.e., C\ + C3 = and C2 + C4 = 0), the continuities 
of -r^r(rw) for m = 0, 1, 2, 3 at r = r , the boundary condition (|73|) and the constraint (|66|). 



Direct calculations reveal that these conditions imply 



1 

1 

a(r o ,0) 
a(r , 1) 
a(r ,2) 
a(r ,3) 








-a(r o ,0) 
-a(r , 1) 
-a(r ,2) 
-a(r ,3) 
a(i? ,0) 
R a(R , 3) - a( J R , 2) 
-^R o a(R o ,l) + ^Si(R o ,0) 










c 2 







c 3 




-pro^o 
-pit 






c 5 







c 6 







c 7 







C 8 








(74) 
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where the 1x4 row vector a(r, m) is given by 



a(r, m) = [&? exp(k + r) , k™ exp(/c_r), (-/c + ) m exp(-/c + r), (-/c_) m exp(-/c„r)]. 



Note that the last row of equation ( |74|) follows from setting the ball Bn to be in equa- 
tion (^) and the constraint (|66|) . Further, from equation (|~|) wc have rcf){r) 
and hence 



la or- 1 a 



f-fgr + C+C^eMk+r) - ex V (-k + r)} 

+C_C2[exp(/i;„r) — exp(— fc„r)] if r < r , 

C+C5 exp(/c+r) + C^Cq exp(A;_r) 

+C + C , 7 exp(— A; + r) + C_C 8 exp(— k_r) if r > r , 



(75) 



where, by equations fl68|) and (0), C± = — Xk 2 z /2a. Therefore, the electrostatic contribution 
to the defect energy is given by 



2 J M 3 2a 



—A f —n^XprQ 
4a 



3a 3 



C+C, 



l7i[k + r cosh(A; + ro) — smh(k + r )} 



k: 



87r[fc_ro cosh(/;;_ro) — sinh(/c_ro)] 

kJ 



• (76) 



We remark that the algebraic equations ( |74|) determine the constants [Ci, • • • , Cg] uniquely 
Analytical expressions of these constants are desirable but impractical to write them down. 
In the limit Rq — > +oo, we find 



Ci 



pl^k 2 _{l + k + r ) exp(-/c + r ) 



2kj r {k 2 _ ■ 

, c 7 



c 9 



-plfa\{l + /c_r ) exp(-/c_r ) 
2fc_(fc+ - k 2 _) 



plQk_[k + r cosh(/c + r ) — sinh(/c + r )] 
M*;2 - fci) 



(77) 



p^o^+[^- r o cosh(/c_r ) — sinh(fc_ro)] 
k_(kl-k 2 _) 



For general cases with finite Rq, which represent computations on a finite simulation cell, 
we resort to numerical solutions. In particular, we are interested in estimating the error 
incurred in the defect energy from using a simulation cell, and its dependence on the cell- 
size. To this end, we have conducted a periodic calculation on a unit cell of FCC lattice for 
aluminum using a real-space formulation for OFDFT and a finite-element discretization of 
the formulation flOavini ct al. , 2007b| ). In our simulation, we used the TFW family of kinetic 
energy functionals with A = | and a modified form of Heine- Abarenkov pseudopotential for 
aluminium flOoodwin et al.| , |1990| ). We subsequently estimate the constants a, /3,j from our 
numerical calculations to be 



a 



0.1629, p = -0.0509, 7 = 0.9449. 



We now estimate the cell-size effects in the electrostatic contribution to the energy of a defect 
that is representative of a vacancy. A reasonable choice for the length scale of a vacancy is 
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r o = a o/2, where ao is the lattice parameter for aluminum which is computed to be 7.5 
a.u. Using equations (ff^)-(ff^), we numerically solve for the electrostatic contribution to 
defect energy. Figure [I] shows our estimate of cell-size effects from finite cell simulations. 
As is evident from these results, Rq = 6r = 3a is necessary for the approximation errors 
from finite cell-size studies to be within 1% of the defect energy — a threshold representative 
of chemical accuracy. In typical electronic structure simulations this Rq corresponds to a 
simulation cell with 6x6x6 FCC unit cells containing 864 aluminum atoms. This estimate 
is in close agreement with recent cell-size studies on vacancy formation energies conducted 
in |Gavini et al.| fl2007a|) , where about 10 3 atoms were required for the cell-size effects in 
defect formation energy to be within O.OleV. We note that despite the exponential decay in 
the electronic fields, cell-size effects are significant, even for a simple defect like vacancy. In 
the more accurate models of density functional theory, like the Kohn-Sham formulation, the 
decay in electronic fields is known to be slower and hence cell-sizes beyond those considered 
in previous electronic structure studies may be needed for an accurate study of defects. 




23456789 

R o /r o 

Figure 1: Cell-size effects showing relative error in the electrostatic contribution to the defect 
energy from finite cell calculations. 

We now consider the elastic contribution of the defect energy, i.e., the first term on 
the r.h.s. of equation (p0|), which is a standard calculation and provide it for the sake of 
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completeness. For simplicity, we assume that the stiffness tensor of the crystal, defined by 
(PI), is isotropic and that the "eigenstress" B is dilatational. Let /i be the shear modulus, k 
be the bulk modulus, and B = (TqI (I is the identity matrix). Based on the Eshelby's solution 
(Eshelby 1957), we find that the displacement is given by 



\Q x r 2 + 9 if |x| < r , 
\® 2 r 2 + 3 ± if r < |x| < Ro, 



where 61,62,63 G IR are constants to be determined. Indeed, by direct calculations we 
verify that the function y given by the above expression satisfies the last of equation ( |51~D 



inside the ball r < r$ and inside the annulus region ro < r < Rq. Across the interface r = ro, 
the continuity of y, the continuity of traction and the boundary condition y = at r = Rq 
imply 

©1 = ©2 - ©sAi pcr + 3k©x = 3kB 2 + 4/i6 3 /r 3 , 6 2 - Q 3 /R 3 = 0. 
Direct calculation reveals that 

pa r 3 Q 
1 Afi + 3K [ R 3 

Therefore, the elastic contribution to the defect energy is given by 

The elastic contribution of the defect energy has a slower asymptotic decay {0{-ks)) in com- 
parison to the electronic contribution and is one other reason to consider large cell-sizes to 
ensure the accurate computation of the energetics of defects. 



7 Extensions 

The form of OFDFT energy we considered for the multiple scale analysis in prior sections 
represents an orbital-free model with TFW kinetic energy functionals without exchange and 
correlation terms. In this section we comment on other general forms of energies that are 
widely used in OFDFT computations. We remark that the multiple scale analysis is indepen- 
dent of the form of the non-linear term f(u) appearing in equation (JJ), and thus including 
the exchange and correlation energies does not affect the analysis or the derived expressions. 
However, the non-local kernel energies can not be represented by a local function of the form 
f(u), and we now present the extension of our analysis to these commonly used kinetic energy 
functional forms. 

The OFDFT formulations employed in numerical studies widely use functional forms for 
kinetic energy that are non-local in real-space, called kernel energies, which are considered 
to be more accurate than the local TFW functionals (cf. equation (fj)). We refer to |VVang 



fe Teter| 01993 ); ISmargiassi fc Madderj ( |l99l) ; Wang et al.| (|T998| , |1999|) for further details 



on these models. We also remark that recent analysis ( [Blanc fc Cances] , [2005 ) has shown 
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that some of the proposed models lack global stability and can pose serious numerical issues. 
For the sake of completeness, we briefly discuss the multiple scale analysis of these non-local 
kernel energies. The functional form of these kernel energies is given by 



E Ker (u) 



J J p{u{x))K(\x-x.'\)q(u{x'))dxdx', (7t 



where p(u),q(u) are functions whose specific form depends on the particular flavor of the 
OFDFT model, and the total energy is given by 

E(cj>, u; b y ) = J [/(«) + \ | Vu| 2 - 1 1 V0| 2 + (u 2 + b y )<f>] rfx + E Ker (u). (79) 

We define the following potentials which will be used to reformulate the non-local kernel 
energy given by equation ( fTlf) into a local variational problem: 

V p (x) = J K(\x-x!\)p(u(x!))dx!, V q (x) = J K(\x-x!\)q(u(x!))dx!. (80) 

Taking the Fourier transform of the above expressions we obtain 

V p (k) = K(\k\)p(k), V q (k) = K(\k\)q(k). (81) 



Following the ideas developed in |Choly fc Kaxirasj ( |2002| ), K can be modeled to good 



accuracy using a sum of partial fractions of the form, 

j = l -vj 

where Pj, Qj, j = 1. . .m are constants which are fitted to best reproduce i^(|k|). These 
constants can possibly be complex, but appear in pairs with complex conjugates. Substituting 
this approximation for K in equation (|8~ID and taking the inverse Fourier transforms, we 
obtain a system of coupled partial differential equations with possibly complex coefficients 
given by 

' -AV PJ + Q,jV m + PjA P (u) = j = 1 . . . m, 
-AV gj + QjV qj + P 3 Aq(u) =0 j = 1 . . . m. 

where V P j and V q j are the inverse Fourier transforms of nriq^g P(k) and . q(k) respec- 

tively for j = l...m. Further, V p w j V pj> V g ~ J2j V qj- B Y defining <p pj = V pj - Pjp(u) 
and <p q j = V q j — Pjqiu) for j = 1 . . . m, equation (|83f) can be rewritten as 

-Aip pj + Qjifpj + PjQjp(u) =0 j = 1 . . . m, 
-Aif qj + QjVqj + PjQjqiu) =0 j — 1 . . . m. 
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The kernel energy, E Ker , can now be expressed in a local form in terms of the potentials 
<fpj, ifgj, or equivalently as a local saddle point problem: 



1 P 1 /* 

E Ker {u) = minmax { V] — — / Vip pj ■ Vip qj dx + — / (p pj (p qj dx 

fpj Vgj < r jLJj J Y J Y 

+ J tp q jp(u)dx + J (p P jq(u)dx. + Pj J p{u)q{u)dx^ . (85) 

We note that variations with respect to tp p j and tp g j return the Euler-Lagrange equations 
given by equation and the saddle point problem correctly represents, within the ap- 

proximation (|82"D, the kernel energy and its functional derivatives. 

We now decompose the potential fields (ip P j,{p q j) into a predictor {<p P j ,<fqj ) and a cor- 
rector (ip p j c ,ipgj c ), and expand the corrector fields using a two-scale expansion given by 

= *) + wk( x > *) + ••■ , (g6) 

= ^>( x ' x ) + W«- C (x,x) + ••• . 

Following on similar lines as in section |, we obtain the following expressions for j = 1 . . . m 
from the leading order terms of the expansion in equation 




Vs4(x, x) = and V*< c (x, x) = 0, (87) 
Vx^ c (x, 5c) = and V*^- (x, x) = 0. (88) 

Thus, the corrector fields in their leading and first order are independent of the fast variable 
representing the lattice length scale. The governing equations for <^° Jc (x) and <£>°- c (x) are 
given by, 



where 



-A<p° pjc + Q jV % c + £ P (F>° = on Y, 
- A<, + Q ]V % c + £ 3 (F>° = on Y, 

£ P (F*) = / p'(n p (F*,x))rfx, f,(F*) = / g'(n p (F*,x))rfx. 



59) 



F*U J F*U 

Further, the governing equations for {(fP c ,u° c ) are given by 

' A0° + 2a(F*)u° + b c = on y 

-A M ° + 27(F*) M ° + 2a(F*)0° + ££i fe(F*)< c + £ ? (F>y =0 on y 



(90) 
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where 



7(F*) = 7 (F*) + - (%*(*") + X qj (F*) + ^i(F*)) : 



Xpj 



;f*) = / i 

J F*U 



X*i(F*) 
^•(F*) 



F*t/ 



/(« p (F*,x))^(F* ) x)dx, 
g"(n p (F*,x))^ p (F*,x)rfx, 



F*t/n 



p"(up(F*, x))g(«p(F*, x)) + 2p'(u p (F*, x))g'K(F*, x)) 



+p( Mp (F*,x))g"( Mp (F*,x)) 



c/x. 



Finally, we comment that the results obtained with OFDFT as the model theory are 
equally valid for the field formulations that describe empirical interatomic potentials pre- 
sented in Iyer fc Gavini ( 2010 ). We note that the field formulation presented in [Iyer &: 
Gavini| ( |2010[ ) result in a system of coupled linear partial differential equations which repre- 



sent a special case of the non-linear governing equations describing OFDFT. 



8 Summary 

The main idea behind the quasi-continuum reduction of field theories is the coarse-graining 
of corrector fields in the formulation using an unstructured finite-element triangulation. In 
this work we have presented a formal mathematical justification that supports such a coarse- 
graining, and places the quasi-continuum reduction of field theories on a firm mathematical 
footing. In particular, we have demonstrated using perturbation method and multiple scale 
analysis that the corrector fields do not exhibit fine-scale (atomic-scale) oscillations in the 
leading order, which allows for the coarse-graining of these fields. Further, we have derived 
the homogenized equations that govern the macroscopic far-field nature of these corrector 
fields, and using Fourier analysis we have estimated their far-field asymptotic behavior. In 
the case of orbital-free density functional theory with TFW kinetic energy functionals, the 
electronic fields comprising of the electrostatic potential and electron density are found to 
exhibit an exponential decay. 

Using the computed asymptotic behavior of these corrector fields, we have estimated the 
errors incurred in the computation of defect energies using finite cell simulations. Although 
the electronic fields exhibit an exponential decay, our analysis shows that cell-sizes of the 
order of 10 3 atoms are required for an accurate computation of defect energies, which is in 
keeping with recent cell-size studies conducted in |Gavini et al.| (|2007a|) . We note that in 



the more accurate versions of density functional theory, like the Kohn-Sham formulation, the 
decay in electronic fields is known to be slower. Further, the asymptotic decay in elastic fields 
is much slower than electronic fields and this effect can become very significant for stronger 
defects like dislocations. This suggests that larger cell-sizes than those that are typically 
used in electronic structure calculations (~ 100 atoms) are needed for an accurate study of 
defects in materials. 
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A priori estimates on the asymptotic behavior of corrector fields from this work can be 
used to determine the optimal coarse-graining rates for finite-element triangulations in the 
quasi-continuum formulation of field theories, and presents itself as a worthwhile future di- 
rection to pursue. Further, developing the quasi-continuum reduction of Kohn-Sham density 
functional theory and an analysis of this formulation is an important open problem, which 
is the focus of our future work. 
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